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^ ■ The origin of stochastic fluctuations in gene expression has received considerable 
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this paper, we propose a stochastic model for gene expression regulation including 
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attention recently. Fluctuations in gene expression are particularly pronounced in 
cellular systems because of the small copy number of species undergoing transitions 
between discrete chemical states and the small size of biological compartments. In 



several binding sites, considering elementary reactions only. The model is used to 
investigate the role of cooperativity on the intrinsic fluctuations of gene expression, 
j> ■ by means of master equation formalism. We found that the Hill coefficient and 

t-h ; 

, the level of noise increases as the interaction energy between activators increases. 

in 

Additionally, we show that the model allows to distinguish between two cooperative 
■ binding mechanisms. 

o\ ' 
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I. INTRODUCTION 

All chemical reactions have intrinsic fluctuations that are inversely proportional to the 
system size. Such fluctuations are particularly pronounced in gene expression. At the tran- 
scriptional level, gene expression is mainly controlled by the cis-regulatory system (CRS) 
and transcription factor (TF) proteins that bind specifically to DNA sites lj. The TFs 
influence the transcription rate by interacting with other transcriptional components (RNA 
polymerase, TATA-binding protein, etc.). Like any molecular interaction, the binding of 
TFs to the regulatory sites is a stochastic event rendering the transition between states of 
the CRS a stochastic process. This source of noise is known as intrinsic noise in gene ex- 
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pression regulation, to distinguish it from that produced by other influences such as random 
fluctuations in nutrients, cell division or regulatory inputs to the transcriptional machinery, 
known as extrinsic noise 

aa. 

There is a broad variety of different CRS motifs that underly such regulation. The diver- 
sity of CRS range from simple ones to more complex motifs that include dozens of regulatory 
sites, some of them organized in clusters or tandems This cluster organization points 
to cooperative effects in the gene regulatory process, because proteins rarely seem to bind 
to DNA without interacting with other DNA-binding proteins. Despite this complexity, the 
bulk of stochastic models for gene regulation are based on transitions between two promoter 
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states (active and inactive) and, recently, more complex models have been explored 
All these models approximate the transcriptional control by using a regulatory expression 
function (Hill function in 6, 7, 8, 9] or an ad-hoc function to fit the model to the experimen- 

fin 

tal data in |4|, |5(). The approximation assumes that changes in the levels of TF are reflected 
instantaneously in the transcription rate. Although this approximation could be reasonable 
to study the static deterministic behavior of transcriptional regulation {si, [9} , it could leads 
to a significant underestimation of transcriptional noise 10j. Consequently, these models 
cannot accurately describe how the overall regulatory process affects noise expression. In 
this article, we propose a theoretical model of transcriptional regulation that considers a 
CRS with several regulatory binding sites for activating proteins. All transition rates be- 
tween CRS states follow the law of mass action for elementary reactions. In this way, our 
model accounts for the fact that the expression response is determined by the dynamics of 
CRS. 



II. THE MODEL 



We are interested in exploring how the molecular interaction affect the cooperativity and 
the fluctuations level of the gene expression. In this sense, we found that stronger interaction 
between activators increases the level of noise expression. In our model the transcriptional 
regulation is assumed to be a stochastic process in which the regulatory system makes 
transitions between different states. The model includes N regulatory binding sites for 
the same TF (Fig. 1 illustrates the case with three regulatory binding sites). The states 
s = l,2,...,A r + l represent, respectively, states with 0,1, ... ,N binding sites occupied by 
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TFs. The states s > N + 2 correspond to the transcriptional complex formation, where all 
components required for transcription are assembled on the CRS. Once the core transcription 
apparatus is formed, the synthesis of one mRNA copy begins. 
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FIG. 1: Kinetic regulatory scheme. All parameters shown in this schematic diagram are constants. 



The working hypothesis is that TFs bound to DNA alter the probability of transcriptional 
complex formation. Consequently, states s < N + 1 are characterized by different rates 
of transcriptional complex formation. For simplicity, we consider that some bindings are 
sequential, i.e., TF does not bind or unbind after transcriptional complex formation and 
transcriptional complex does not assemble before TFs bound to DNA site. Additionally, we 
consider that the sites are functionally identic, i.e., the model does not distinguish among 
states with the same number of TFs bound to regulatory binding sites. Thus, in our model, 
the states of CRS are related more to the occupancy number rather than to the binding 
status of each site. This additional simplification reduces the number of states accessible 
to the CRS and allows us to explore the role of cooperativity on noise expression without 
considering a combinatorial number of states. The model assumes that mRNA is synthesized 
at a rate which depend on the state s. mRNA is considered to degrade linearly with rates 
7- 



As other authors , we used the master equation approach to derive the average 

behavior of mRNA level, as well as its fluctuations. For this system, the state is specified 
by two stochastic variables: the state of the CRS s, and the number of transcripts m. We 
can write the probability to find, at any given time t, the system in the state (s, m) as a 
vector P m (t) = (P^ (t) , (t) , . . . , P m N+1 (t)j. The time evolution for this probability is 
governed by the following master equation: 

2N+1 

P s = V t P r + a (P s , - P s ) 

r=l 

+ 7 [(m+l)Ci-<] C 1 ) 

where t s ^ r are the elements of the transition matrix T and represent the transition probability 
per time unit from state r to state s, a s are the elements of a vector and correspond to 
the transcription rate of state s. The first term on the right hand side of the equation ([1]) 
describes the CRS dynamics, while the second and third terms correspond to the production 
and degradation of mRNA, respectively. 

We are mainly interested on both the gene expression level and its fluctuations. The first 
is measured through the first moment of the number of transcripts m, and the last through 
the corresponding variance, related to the second moment: 

<™>=E™^ * m = (m 2 )-(m)\ (2) 

m,r 

where the summation limits were suppressed for the sake of readability. We want to remark 
that from now on, every sum over transcript numbers will be from m = to m = oo, while 
the sum over CRS states will be from r = 1 to r = 2N + 1. 

Moments of j-th order can be written in term of its associated partial moments 

(m- 7 ) = ^(m J ) r where (m J ) s = m J P^, (3) 

r m 

Note that the zero partial moments are the marginal probabilities for the CRS to be in state 
s, regardless the number of transcripts, i.e., (m°) s = P s = J2 m P m - 

From Eqs. (CD) we can derive a set of ordinary differential equations for the time evolution 
of the partial moments for any j. As there is no feedback, the system is linear and can be 
solved analitycally. Thus, the time evolution of partial moments for j = 0, 1, and 2 are given 
by 

3 = P S = J2rts,rP r (4) 
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3 = 1 (m) s = ErM m >r + ®sP s ~ l{m) s (5) 
j = 2 ( m ' 2 ) s = Er V(^ 2 )r " 27(™ 2 > s + 2a,(m), + 7 <m) s + a s P s , (6) 

From these partial moments, we can readily find first order differential equations governing 
the time evolution of the mean and variance, 

( m ) = _ 7 ( m ) + ]Ta r P r (7) 

r 

a 2 m = -27 (m 2 ) + 7 (m) + 27 (m) 2 

+ $> r [2(m) r + (1 - 2(m»P r ] . (8) 

r 

From Eq. the steady-state solutions for the mean value of m is 

(m*> = -£« r P;, (9) 

7 ; 

where * denotes the steady state solution. The steady-state solution of the probability vector 
P* corresponds to the normalized eigenvector associated to the zero eigenvalue of the CRS 
transition matrix TP* = 0. 

The steady state solution for the variance follows from Eq. ([8]) 

o 2 ^ = (m*) - (m*) 2 + -£a r (m*) r , (10) 

7 r 

where (m*) r is determined as the solution of the linear equation 

(V ~ T^s,r) (m*) r = -a s P*, (11) 

r 

where <5 Sjr is the Kronecker delta. The expressions (191- fTTj) are general, in the sense that they 
are valid for iV binding sites in the CRS. In this paper, we have limited the study to the 
case of Fig. 1, i.e., a CRS with N = 3, and a s = a for s > 4 and zero otherwise. We are 
motivated to set N = 3 because the cooperative effects are more apparent for greater N. 
However, an approximation used in the next Section, could not be adequate for higher N. 
The TFs can bind to regulatory sites with a probability proportional to TF concentration 
c (cfci2,c/c23 and C/C34, where kij are the kinetic rates), following the law of mass action for 
elementary reactions. TF unbinding events depend only on the kinetic constants (^21,^32 
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and /C43). In this case the transition matrix T can be written as 
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(12) 



and the associated steady state solutions of the partial probabilities P" involved on Eq. (Q 
are 

cK 2 K 5 



* 

p6 
Pi 



where 



~ 1 + c (K 2 (K 5 + 1)) + C 2 (K 2 K 3 (K 6 + 1)) + C 3 {K 2 K 3 K 4 (K 7 + 1)) 

c 2 K 2 K 3 K 6 

" 1 + c (K 2 (K 5 + 1)) + c 2 (K 2 K 3 (K 6 + 1)) + c 3 (K 2 K 3 K, (K 7 + 1)) 

c 3 K 2 K 3 K A K 7 

= 1 + c (K 2 (K 5 + 1)) + c 2 (i^ 2 K 3 (K 6 + 1)) + c 3 (K 2 K 3 K, (K 7 + 1)) 

fVS — 1,S 



(13) 



for s = 2, 3, 4 and i^, 



for s = 5, 6, 7. Replacing these expressions 



fcs,s — 1 ^ fcs,s — 3 

on Eq. ([9]), we find the explicit form for the steady state expression level of transcripts 
, a K 2 K 5 c + K 2 K 3 K 6 c 2 + K 2 K 3 K 4 K 7 c 3 

TYl —— 

7 1 + K 2 (l + iT 5 )c + ^2^3(1 + K 6 )c 2 + K 2 K 3 K 4 {1 + i^c 3 

Unfortunately, though a closed expression for the variance is obtained it is too long to report 
here. From this point the * will be suppressed from the steady states expressions. 



(14) 



III. BINDING COOPERATIVE MECHANISMS 



As the regulatory sites are assumed to be functionally identic, we can introduce a rela- 
tionship between the TF binding/unbinding when there is no interaction between the TFs. 
Thus, if the probability per time unit that a single TF molecule binds to a regulatory site 
is p, we have k° 2 = 3p, k 23 = 2p, k 3A = p, where indicates that there is no interaction 
between TFs. Similarly, unbinding rates are given by k 2l = q, k 32 = 2q, k% 3 = 3q, where q 
is the probability per time unit that a single TF molecule unbinds from an occupied site. 
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Now, we will assume that the probability for a TF molecule to bind to a given regulatory 
site arises from: i) the free energy of binding a TF to the specific site AGdna ; n ) the 
free energy of interaction between TF molecules bound to adjacent sites AGi. From these 



assumptions and the principle of detailed balance [12j], we are able to find a relationship 
between kinetic rates with and without interactions among TFs. In the case AGi = 0, we 
have 



AGr 



ckl s+1 /k° +hs = e-^~ for a = 1, 2 and 3, (15) 

where h° s s+1 represents the transition rate from state s to state s + 1 when there is no 
interaction between TFs (k° +1 s represents the rate for reverse transition). In general, the 
TF molecules interact with each other, i.e., AGi < 0. Consequently, we have 

^£±i =£ a,^£±l (16) 

_ AGj 

where e = e rt represents the intensity of the interaction between TFs. a s represents 
the number of interactions. We assume that all bounded TFs interact with the new one, 
regardless of their position on the CRS, thus we have a s = s — 1. This means for example 
that in the state s = 4, the third TF interacts with the two bounded TFs. As the cooperative 
effects are more apparent for greater N, we are motivated to use a high number of binding 
sites. However, the assumption a s = s — 1 could be inadequate for s > 4 when N > 3 
due to the greater separation between faraway sites. For this reason, we set N = 3 for 
further calculations. The relationship ffl6l) leaves an extra degree of freedom, because the 
interaction between TFs can increases the binding rate k StS+ i, increasing the ability for new 
TF recruitment for DNA binding, or it can diminishes the unbinding rate k s+ i jS , increasing 
the stability of the TF-DNA bound. The first case will be denoted here as the recruitment 
mechanism (RM), while the second case will be denoted as stabilization mechanism (SM). 
These two mechanisms are not mutually excluding, and they could be acting simultaneously 
in real life, but in order to study their effect on the regulatory response and its associated 
fluctuations, we will consider the alternative cooperativity binding mechanisms separately. 
Thus, using relations (Tl6|) and the relations for binding/unbinding rates, we obtain 

k s>s+1 = e^(4-s)p 

k s+1>s = sq, (17) 
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for the first mechanism, while for the second mechanism we have 



k S)S+1 = (4 - s)p, 
ks+i,s 



(18) 



Table I summarizes the parameter values used in this work (the time unit is min and concen- 
tration is an arbitrary unit.) Binding and unbinding parameters were obtained considering 
p = 0.25, q = 0.75 and e = 6 (which is equivalent to AGi ~ —1.0 Kcal). Cases RM and SM 
have the same TF interaction intensity, while e = when there is no interaction between 
TFs. 



TF binding/unbinding rates 
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SM case 
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0.75 


0.75 


k25 


0.50 
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0.50 
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^36 
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&32 


1.50 
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1.50 


^63 


0.50 


&34 


9.00 


0.25 


0.25 


k i7 


1.50 


h3 


2.25 


0.0625 


2.25 


k 1A 


0.50 




Production and degradation rates 


mRNA (layer III) 


a 


1.50 


7 


0.03 



TC formation rates 



TABLE I: Kinetic parameters values. RM case: recruitment mechanism; SM case: stabilization 
mechanism; eo case: there is no interaction between TFs. The time unit is min and the concentra- 
tion is an arbitrary unit. 



IV. RESULTS 

From the Eq. (fT4j) for the mean, Eqs. (11 OH 111) for the variance, and using the parameters 
values of Table I, we study the response of an activator switch. We consider that the kinetic 
rates for transcriptional complex formation increase linearly with the occupancy number, 
i.e. k SjS+ 3 oc s. With this condition we assume that there is no synergism between TF 



and transcriptional complex formation 



131 ] . This synergism can contribute to the effective 
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cooperativity (data not shown). Figure 2A depicts the average number of mRNA copies (m) 
and Figure 2B depicts the standard deviation (J clS 8b function of the activator concentration 
c, obtained analytically for the three cases: RM case, the interaction between TFs increases 
the binding rates, (solid line); SM case, the interaction between TFs decreases the unbinding 
rates; and (eq) case, where there is no interaction between TFs (dotted line). The mean 
response of RM and SM cases are exactly the same. The regulatory function for both 




c 

FIG. 2: Average number of mRNA copies as a function of c in steady state for three different cases: 
RM (solid line) , SM (overlapping the former case dashed line) , and sq case (dotted line) . See Table 
I for parameter values. (B) Associated standard deviations of the transcript number as a function 
of c (C) Associated noise versus the transcriptional efficiency. 

examples with cooperative binding fit the Hill function / (c) = V ma x/ [1 + { C /Kd) nh ], where 
c is the binding protein concentration, V max is the maximum number of transcripts, is the 
dissociation constant, and is known as the Hill coefficient which determines the steepness 
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of the regulatory function. Both RM and SM cases the Hill coefficient is 1.86, while the 6$ 
case is associated to a response without cooperativity = 1.00. In all cases the fluctuation 
level estimated by the standard deviation a m (see Figure 2B) has a peak centered around the 
dissociation constant Ka. Relative fluctuations are characterized by the normalized standard 

n 

deviation. Analyzing an un- normalized measure can lead to artefactual results [14j. Fig. 
2C depicts the noise rj (defined as a m j (m)) as a function of the transcriptional efficiency m r 
(defined as the ratio between transcription and maximum transcription (m) / (m) max , where 
(m) max is the maximum of (m)) for each case in Figure 2A. 

Figure 2 shows that cooperativity has an effect not only in controlling the expression 
response increasing rih, but also in increasing the relative size of fluctuations. We note also 
that though the regulatory function in RM and SM cases is the same, the mechanism of 
increasing the TF-DNA complex stabilization (SM) is associated to higher level of noise 
(dashed line) than the mechanism involving an improvement in the recruitment ability of 
new TF to the DNA (RM) (solid line). Figures 2 suggest that the two cooperative binding 
mechanisms considered here affect the fluctuation level in a differential way but not the 
regulatory function. Even though this function is altered by the FT interactions, it is not 
possible to distinguish between the alternative mechanisms through the regulatory function 
only. 

We also analyze how the binding and unbinding rates, p and q respectively, affect the 
regulatory function and the fluctuation level. Figure 3A depicts the behavior of the disso- 
ciation constant Kd, as function of the unbinding rate q (left) and as function of binding 
rates p (right), keeping the other rates constant. Kd does not depend on which cooperative 
binding mechanism is acting (the curves are completely overlapped). However, K& increases 
linearly with binding rate q, and it is inversely proportional to the binding rate p. Figure 
3B depicts the behavior of the Hill coefficient n^, as function of the binding/unbinding rate 
(right/left), keeping the other rates constant. We can observe that nh is almost independent 
on these rates. From Fig. 3C, we note that the level of noise is sensitive to the type of 
cooperative binding mechanisms which is. (T max decreases with the unbinding rate more 
slowly in the SM (dashed line) than in the RM (solid line). The difference between the 
two mechanisms diminishes when the unbinding rate decreases, while the maximum value 
of dispersion is not affected when the unbinding rate p varies. Due to the complexity of the 
problem are not able to provide a quantitative proof that SM leads to larger noise than the 
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FIG. 3: Half maximal value of the mean response as a function of the unbinding rate q, and 
as function of the binding rate p. (B) Hill coefficient rih as a function of the unbinding rate q and 
as function of the binding rate p. (C) Maximum value that reaches the standard deviation vs the 
unbinding rate q, and as function of the binding rate p. Left: q varies keeping p = 0.25, Right: p 
varies keeping q = 0.75. For RM mechanism (solid line) and SM mechanism (dashed line). 

RM in all conditions. However, from figure 3C we can see that noise level does not depend 
on the binding rate p, but it depends on the the unbinding rate q. As the interaction energy 
between TFs decreases the unbinding rate in the SM, it is expected that SM has associated 
a higher level of fluctuation than the RM. 

Finally, we show how the interactions between TFs alters both and nh parameters 
of the regulatory function and also the fluctuation level. Figure 4A illustrates how the pa- 
rameters Kd and nh are affected by the interaction intensity e. The Hill coefficient (filled 
circle symbols), scaled on the right vertical axis, increases with e, suggesting that the steep- 
ness of the regulatory function depends linearly on the free energy AG\. Furthermore, the 
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dissociation constant (open square symbols), scaled on the left vertical axis of Figure 
4A, decreases with the interaction intensity. Finally, we found that the fluctuation level 
increases with the interaction intensity. Figure 4B depicts the maximum value of the stan- 
dard deviation a max as a function of e. We have observe that the RM (solid line) is less 
sensitive to the interaction intensity than the SM (dashed line). We want to remark that 
the differences between the recruitment and stabilization mechanisms vanish when there is 
not interaction energy between TF (AG/ = 0). These observation suggest that our model 
predicts that interactions between TFs improve the response of the regulatory system in the 
sense of specificity (higher rih) and sensitivity (lower K^). But, in contrast, the system loses 
accuracy because the noise increases with the intensity of the interaction. 



V. CONCLUSION 



We have shown that a model which includes several binding sites is able to address the 
question of cooperative binding effects on fluctuations. The first moment of m is the same as 
that obtained from thermodynamic models, which depends solely on equilibrium constants. 
Nevertheless, second moments have allowed us to introduce new quantitative insights on the 
TF cooperative binding effects in the cell-to-cell variability. We found that two different 
cooperative binding mechanisms can be distinguished: the RM which increases the ability 
for new TF recruitment, and the SM which increases the stability of the TF-DNA bound. 
In both mechanisms, the Hill coefficient and level of noise increase as the interaction energy 
between activators increases. Only a few kilocalories of binding energy between TFs have 
a dramatic effect on the noise level, which also depends on the acting cooperative binding 
mechanisms. The other hand, the mechanism that reduces the unbinding rates is associated 
to a greater level of noise which is in agreement with two state model 15]. 

Both mechanisms reported here are derived from the thermodynamics relationship used 
in the modeling. This cannot be done in simpler models that use regulatory expression 
function rather than TFs thats bound to several binding sites on DNA following the law of 
mass action. These different mechanisms have not been reported previously. Although the 
proposed model is more complicated than previous, it can also be solved analytically. Thus, 
the model constitutes an adequate frame to discuss the impact of the diverse cooperativity 
mechanisms on the gene expression fluctuations. However, we want to remark that the 
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FIG. 4: (A) Half maximal value of the mean response as a function of e (scale on left axis). 
Hill coefficient as a function of e (scale on right axis). The curves associated to RM and SM 
are completely overlapped. (B) Maximum value that reaches the standard deviation cr max as a 
function of e. Solid line corresponds to the recruitment mechanism, while dashed line corresponds 
to the stabilization mechanism. 
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presented model is limited to intrinsic contribution of noise, i.e. it does not regard the 
fluctuation on the TF concentration and other extrinsic source of noise, which certainly 
contribute to the total noise. 
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